clear
close all
clc

%% Parametros de ajuste de la simulacion
M = 16;  % cantidad de simbolos
k = log2(M); % bits por simbolo
SNRbdB = 0:10; % SNR por bit en dB
nsimbs = 2e6; % cantidad de simbolos que hay para cada valor de SNRbdB en el Montecarlo


%% Conversiones e inicializaciones varias
SNRbveces = 10.^(SNRbdB/10); % SNR por bit en veces
SNRdB = SNRbdB + 10 * log10(k);  % SNR en dB
SNR = 10.^(SNRdB/10);    % SNR en veces
Eavg = 1;
Pes = zeros(size(SNRbdB));  % probabilidad de error por simbolo simulada
Peb = zeros(size(SNRbdB));  % probabilidad de error de bit simulada

%% P(error) teorica:
% Probabilidad de error de bit teorica de una M-QAM
% Uso un puntero a una funcion cuyo parametro es SNR (en veces)
Pet = @(SNR) 4 * (1-1/sqrt(M)) * qfunc(sqrt(3 * SNR/(M-1)));

% Simulacion Montecarlo QAM


modQAM = modem.qammod('M', M, 'SymbolOrder', 'Gray','InputType', 'Integer')
             
demodQAM = modem.qamdemod('M', M, 'SymbolOrder', 'Gray')

const = modQAM.modulate([0:M-1]); % sirve para escalar el resultado luego

for n0 = 1:length(SNRbdB)
    
    %% Primero trabajo con simbolos y despues analizo por bit
    
    sigma = sqrt(Eavg/SNR(n0)); % desvio de ruido
    simbtx = randi(M,nsimbs,1)-1;   % simbolos a transmitir (en 1:M);
    
    scale = modnorm(const, 'avpow', Eavg); % calculo un factor de escala
    
    % la funcion "step" llama al objeto modulador o demodulador
    simbmod = modulate(modQAM, simbtx)*scale; % simbolos modulados para su transmision "OJO, son COMPLEJOS"
      
    % Sumo el ruido (OJO sumo el ruido como si la constelacion tuviera parte compleja)
    simbmod = simbmod + (sigma/sqrt(2)) * (randn(nsimbs,1)+1i * randn(nsimbs,1));
    
    % Demodulo
    simbrx= demodulate(demodQAM,simbmod/scale); % simbolos recibidos
    
    % Probabilidad de error de simbolo simulada
    Pes(n0) = nnz(simbtx~=simbrx)/nsimbs;
    
    %% Ahora analizo los errores de bit (el Gray lo hizo el modulador)
    
    bitstx = de2bi(simbtx); % convierto los simbolos transmitidos a bits
    bitsrx = de2bi(simbrx); % convierto los simbolos recibidos a bits
    % comparo los bits y cuento cuantos difieren
    Peb(n0) = nnz(xor(bitstx,bitsrx))/k/nsimbs; 
    
end

%% Graficos
% Probabilidad de error de simbolo
figure(1)
semilogy(SNRbdB,Pet(SNR),'-dr','linewidth',1.5);
hold on
semilogy(SNRbdB,Pes,'-o','linewidth',1.5);
grid on
box on
xlabel('SNR por bit (dB)');
ylabel('P_e')
legend('Curva teorica','Montecarlo')

% Probabilidad de error por bit
figure(2)
semilogy(SNRbdB,Pet(SNR)/k,'-dr');
hold on
semilogy(SNRbdB,Peb,'-o','linewidth',1.5);
grid on
box on
xlabel('SNR por bit (dB)');
ylabel('P_b')
legend('Curva teorica','Montecarlo')
